DaTa Preparation NBS

A framework to assess the potential ecosystem service gains from restoration across five landscapes.

Purpose

Estimate the potential value of ecosystem services provision gain if the target surface is distributed randomly across the intervention areas.

The ecosystem data included here is derived from Chaplin-Kramer et. al (2022) Mapping the planet’s critical natural assets

1. Coastal Protection. Unitless measure, refers to a vulnerability index. InVEST Coastal Vulnerability Model

2. Nitrogen Export. Derived from the Nitrogen retention modeled using the InVEST Nutrient Delivery Ratio. Expressed in kg/pixel/year

3. Sediment Retention. Derived using InVEST SDR: Sediment Delivery Ration. Values in ton/pixel/year

4. Pollination. Derived from InVEST SDR: Pollinator Abundance Model. Units expressed on equivalent people feed. More information in Chaplin-Kramer, et al. 2022

5. Nature Access represented as the number of people within 1 hour travel of natural and semi-natural lands (Chaplin-Kramer et al, 2022).

   Goals:

• Identify High-Value Pixels: Focus restoration on areas where ecosystem service gains are maximized.

• Estimate Total Gains: If restoration is distributed randomly over a target area, what is the expected ecosystem service value?

• Represent uncertainty through confidence intervals and spatial distributions of selected pixels.

Methods:

Assumptions

• All ecosystem services are equally valued.
• Value is a constant function of service quantity.
• Uniform likelihood of restoration across the intervention area unless constrained by inputs (e.g., masks, distance).
• Independence of pixel values (no spatial autocorrelation considered in sampling)
• Sampling provides a robust method for estimating random distribution effects.

Approach:

• Normalize ecosystem service rasters to a 0-1 scale, add them, and identify high-value areas.
• Randomly sample pixels to estimate service values under random distribution scenarios, iterating multiple times for robust confidence intervals. -Address spatial distribution issues (e.g., clustering) later using additional weighting factors (e.g., downstream beneficiaries).

Combined ES raster production:

Preliminar first step to deal with raster data sets representing different variables.

Limitations: - Implicitly assumes equal importance of ecosystem services. - May overemphasize areas with one dominant service, especially pollination, which is contingent to agricultural production. - Fails to account for geographically weighted factors, plan to address.

Random Sampling:

Suited methodology when no location-based prioritization is assumed. Repeated sampling process assuming Central Limit Theorem. Summarize results with total estimated ES gains and confidence intervals.

1. Prepare Environemnt

Load libraries, functions and prepare template data

1.2 Load Templates

This loads and prepares template data used to align

# add backgrounds/templates to align 
path_lc <- here('ESA_LC') 
# load reclassified land cover map
tf <- file.path(path_lc, list.files(path_lc, pattern= "rec"))
lc <- lapply(tf,rast)
# create rcl matrix
rcl <- matrix(c(
  0, Inf, 0   # Any value greater than 0 becomes 1
), ncol = 3, byrow = TRUE)

#create background pixels, subsitute all by 0
tmp <- lapply(lc, function(r){
  r <- classify(r[[1]], rcl)
})

rm(lc)

2 Extract the data

2.1 List the Restoration Layers

#Restoration:
tiffes <- file.path(here("cropped_raster_data" ), list.files(paste0(here("cropped_raster_data")),pattern= '.tif$'))
#Remove points, risk reduction baselin    
tiffes <- tiffes[c(1:5,51:55,71:75,81:90,106:110)]
#tiffes <- tiffes[-c(6:10)]
#tiffes <- tiffes[c(6:10)]

2.2 Normalize and combine data layers

The combined raster is a proxy for the total ES provided(with some assumptions that need to be eventually reviewed and refined).

2.3 Load Adjusted Griscom Restoration Data and apply mask

Use the adjusted Griscom restoration data (deals with pixel values encoding issues in the original data, had to get the metadadata for that). We apply this mask to isolate only the pixels that have been identified as Restoration Potential in Griscom’s Dataset on Global priority areas for ecosystem restoration.

Intervention Areas & Combiend Restoration Potential

3 Calculate Restoration Potential output.

tiffes <- file.path(here("restoration_combined"), list.files(paste0(here("restoration_combined")),pattern= 'ES_msk'))
baseES_m <- lapply(tiffes, rast)
tiffes

[1] "/home/jeronimo/OneDrive/WWF_nbs_op/restoration_combined/BRAZIL_ES_msk.tif"    
[2] "/home/jeronimo/OneDrive/WWF_nbs_op/restoration_combined/MADAGASCAR_ES_msk.tif"
[3] "/home/jeronimo/OneDrive/WWF_nbs_op/restoration_combined/MEXICO_ES_msk.tif"    
[4] "/home/jeronimo/OneDrive/WWF_nbs_op/restoration_combined/PERU_ES_msk.tif"      
[5] "/home/jeronimo/OneDrive/WWF_nbs_op/restoration_combined/VIETNAM_ES_msk.tif"   

3.1 Brazil Espiritu Santo

3.1.1 Prepare Template Brazil Espiritu Santo

This just create a background of Zeroes to use as template for reprojecting (this is necessary because the access layer whic is on a different crs. Not ideal, but need to move forward)

3.2 Madagascar (Not Necessary To Mask for te Interventions)

3.3 Yucatan

3.3.1 Prepare Template Yucatan

This just create a background of Zeroes to use as template for reprojecting (this is necessary because the access layer whic is on a different crs. Not ideal, but need to move forward)

3.3 Peru

3.4 Vietnam

Pending…

4 Sampling and Extracting values

This will be dealt on a two part basis: 1.Identify the hectares that yield the maximum aggregated restoration values withing the Griscom restoration potential pixels. Contingent on the entry assumptions 2. Model the estimated values of ES gains for the services based on a random sampling. This second part has implicit a couple of assumptions:

• We are not considering any spatial configuration (landscape metrics) aspects affecting the total value
• Eventually, the model can be refined by incorporating additional parameters (e.g minimum distance to boundaries, distance between points, topography or inhabited areas (some of it is implicit in the input data, but needs to be confirmed).
• Coastal risk protection only occurs at the coast, any random sampling performed will haveto consider this. This can be adjusted using weights.

This approach utilizes a stratified random sampling technique to estimate the mean values of different bands in a multi-band raster dataset. Stratified random sampling is a probability sampling method where the population is divided into homogeneous subgroups called strata, and random samples are taken from each stratum. In this case, the strata are defined by the spatial extent of the raster dataset, and the pixels within the raster represent the individual sampling units.

Methodology

Define the Area of Interest: The first step is to define the area of interest (AOI) within the raster dataset. This AOI represents the spatial extent from which the samples will be drawn. In this specific case, the AOI is defined as 30,000 hectares.

Calculate the Sample Size: Based on the resolution of the raster and the desired AOI, the required number of pixels to be sampled is calculated. This ensures that the total area covered by the sampled pixels corresponds to the defined AOI.

Perform Stratified Random Sampling: The spatSample() function from the terra package in R is used to perform stratified random sampling. This function allows for random sampling of pixels within the defined AOI, while also excluding pixels with “NA” values in all bands, ensuring that only valid data points are included in the analysis.

Repeat Sampling: To reduce potential sampling bias and improve the accuracy of the estimates, the random sampling process is repeated multiple times (in this case, 30 times). This is analogous to the concept of bootstrapping, where repeated sampling with replacement is used to estimate the sampling distribution of a statistic.

Calculate Summary Statistics: For each repetition, the mean, standard deviation, and 95% confidence intervals are calculated for each band in the raster dataset. This provides a measure of the central tendency and variability of the sampled data.

Synthesize Results: The results from all repetitions are combined to calculate an overall mean and confidence interval for each band. This provides a more robust estimate of the expected values, effectively correcting for potential outliers and reducing the influence of individual sample variations.

Theoretical Justification

The use of stratified random sampling is justified as it ensures that the sample is representative of the entire population (i.e., all pixels within the AOI). By dividing the population into strata and sampling from each stratum, this method reduces the variability of the estimates compared to simple random sampling. This is particularly important when dealing with spatial data, where there may be inherent spatial autocorrelation or heterogeneity.

The repeated sampling approach further enhances the robustness of the estimates by providing a distribution of possible values. This allows for the calculation of confidence intervals, which provide a measure of the uncertainty associated with the estimates. By synthesizing the results from multiple repetitions, the overall estimates are less susceptible to the influence of individual sample variations and provide a more accurate representation of the true population values.

References

Cochran, W. G. (1977). Sampling techniques (3rd ed.). John Wiley & Sons. Lohr, S. L. (2010). Sampling: Design and analysis (2nd ed.). Brooks/Cole. 1
1. Bivand, R. S., Pebesma, E. J., & Gómez-Rubio, V. (2013). Applied spatial data analysis with R (2nd ed.). Springer.

4.1 Espirito Santo

4.1.1 Prepare Data Esp Santo.

Align and resampe all the raster datasets to the same crs and origin.

4.1.2 Run Sampling and Syntetize results - Espirito Santo

Here, we are sampling the data for obtain estimate of the expected restoration gains assuming randomly selected pixels extracted form the potential restoration areas. Again, some assumptions will have to be reviewed, but this is an initial assessment. The targeted intervention area is 30.200 ha.

4.1.3 Plot Results Espirito Santo

4.2 Yucatan

The approach is the same as in Brazil. Two main intervention areas were used tro build the final polygon: the jaguar corridor and ### 4.2.1 Prepare Data Yucatan. Align and resample all the raster bands to the ame crs and origin.

4.1.2 Run Sampling and Synthetize results - Yucatan

4.1.3 . Plot Results Yucatan

Next Steps:

Run the Analysis for Madre de Dios, incorpporate beneficiaries/population related data, as the target metric is not the sum of the values.